#!/usr/bin/perl

my $option_file = "panam.ini";
chomp($option_file);

die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";

my $usearch;
my $NGS_id_Results;
my $panam_output = "panam_output";
# my $parsing;
my $query_seq; 
my $user_file;
my $preprocess_output = "preprocess_output";
my $seq_F;
my $seq_R;
my $primF;
my $primR;
my $trim =0;

my $Reference = "Reference";



####
my $dom;
my %dom;
my $aff_norm; 
## VERSION 3.02 ##
my $nbrBar; my $ficBar;
my $barIn; my $barcode;
#### V3.02

while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}
	
# get input options line
        else {
		
		if($line=~m/454_RUN_IDENTIFIER/gi) { 
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; 
			}
	 		die "\n\tOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( $NGS_id_Results ne "");
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "\n\tCould not create Directory $NGS_id_Results $!\n"; }
		}

# 		if (-d $NGS_id_Results."/".$panam_output) { `rm -r $NGS_id_Results/$panam_output` } 
		if (!(-d $NGS_id_Results."/".$panam_output)) { mkdir $NGS_id_Results."/".$panam_output || die "\n\tCould not create Directory $NGS_id_Results/$panam_output $!\n"; }


		if($line=~m/QUERY_SEQUENCES/gi) {
			$line=~s/QUERY_SEQUENCES//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi; 
			if($line) {
				$query_seq = $line; 
			}
			chomp ($query_seq);
	
		}

		if($line=~m/BARCODE_FILE/gi) { 
			$line=~s/BARCODE_FILE//gi; ;
			$line=~s/\t//gi; 
			$line=~s/^\s+//gi; 
			$line=~s/\s+$//gi; 
			if($line) {
				$barIn = $line; #print "$barIn\n";
			}

	
			if (!(defined $barIn ) or ($barIn eq "")) {
				$barcode =0;
				$nbrBar = 0;
				########ajout 13/5
				$ficBar = "no barcode file" ;
				#############fin ajout 13/5
			}
			elsif ((defined $barIn) and ($barIn ne "")){
				$barcode = 1;
				##########" ajout 13/5
				$ficBar = $barIn;
				########fin ajout 13/5     

				# See if the barcode file exist and contains something
				die "\n\n\tCannot find barcode file: $barIn. Check panam.ini.\n\n" unless (-e $barIn);
				die "\n\n\tBarcode file, $barIn, appears to be empty!\n\n" if (-z $barIn);

				open (BB, $barIn);
				my @barIN= <BB>;
				close BB;
				foreach my $e (@barIN) {
					if ($e !~ /\t/) {die "\n\n\tYour bar code file $barIn does not seem to be in the right format.\n\n";}
					my @ty = split(/\t/, $e);
					if (($ty[0] =~ /^\d/) or ($ty[0]=~ /_/) or ($ty[0]=~ /\//)) { die "\n\n\tThe bar code Ids in $barIn should start with no number and contain no underscorses.\n\n";}
					if ($ty[1] !~ /[ATCGU]/i) {die "\n\n\tYour bar code file $barIn does not seem to be in the right format.\n\n";}
				}
				############### fin modif + ajout

				my $nbr = `wc -l $barIn`;
				if($nbr =~ /(\d.?)\s/) { $nbrBar = $1} #### modif 25/5 w --> d 
			}
		}
		#### V3.02

		if($line=~m/USER_FILE/gi) {
			$line=~s/USER_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$user_file = $line; 
			}
			die "\n\n\t User file value missed. Check panam.ini.\n\n" unless (($user_file ne ""));
			die "\n\n\t Invalid user file value!\n\n" unless (($user_file eq "own") or ($user_file eq "pooledSamples") or ($user_file eq "eachSample"));	
		}

		if($line=~m/FORWARD_PRIMER_NAME/gi) {
			$line=~s/FORWARD_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primF = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_NAME/gi) {
			$line=~s/REVERSE_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primR = $line; 
			}
		}

		if($line=~m/FORWARD_PRIMER_SEQUENCE/gi) {
			$line=~s/FORWARD_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_F = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_SEQUENCE/gi) {
			$line=~s/REVERSE_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_R = $line;
			}
		}
			
		############### normalisation
		if($line=~m/NBR_SEQ_NORM/gi) {
			$line=~s/NBR_SEQ_NORM//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$aff_norm = $line; #print "aff --- $aff_norm\n";
			} 
		}
		
		#####################################

		if($line=~m/DOMAIN/gi) {
			$line=~s/DOMAIN//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) { 
				$dom = $line;
				if (defined $dom ) { #push (@dom, $dom);
					$dom{$dom} = 1
				}
			}
			foreach my $d (keys %dom) {
				die "\n\n\t Domain name is not correct. Check panam.ini.\n\n" unless (($d eq "bacteria") or ($d eq "eukaryota") or ($d eq "archaea"));
			}	
		}
 	}
}


die "\n\n\t Invalid number of sequences to pick for the normalization! Check panam.ini.\n\n" if ((defined $aff_norm ) and ($aff_norm !~ /^\d+$/));
if (!(defined $aff_norm)) {$aff_norm = '200' }


## VERSION 3.02 ##
##################### Ajout JCC avril 2013 : by sample order

my $socat = `cat -n $barIn | awk -F"\t" '{ print\$1"\t"\$2 }'`;
my @socat = split("\n", $socat);
my %so = ();

foreach my $line (@socat) {
    my ($val, $key) = split("\t", $line);
	
    $so{$key} = $val;
}

#########################################################
#### V3.02




########################### ajout position

my $forw_regexp; my $reve_regexp;
my @tab_forw=split(//,$seq_F);
$forw_regexp=regexp(@tab_forw);

my $rev=reverse($seq_R);
my @tab_reve=split(//,$rev);
$reve_regexp=regexp(@tab_reve);
$reve_regexp=~tr/ACGURMBD/UGCAYKVH/;

sub regexp{
	my $regexp;
	foreach (@_){
		if ($_ eq "T"){
			$regexp.="U";
		}
		elsif ($_ eq "R"){
			$regexp.="[AGR]";
		}
		elsif ($_ eq "Y"){
			$regexp.="[CUY]";
		}
		elsif ($_ eq "M"){
			$regexp.="[CAM]";
		}
		elsif ($_ eq "K"){
			$regexp.="[UGK]";
		}
		elsif ($_ eq "W"){
			$regexp.="[UAW]";
		}
		elsif ($_ eq "S"){
			$regexp.="[CGS]";
		}
		elsif ($_ eq "B"){
			$regexp.="[CUGB]";
		}
		elsif ($_ eq "D"){
			$regexp.="[AUGD]";
		}
		elsif ($_ eq "H"){
			$regexp.="[AUCH]";
		}
		elsif ($_ eq "V"){
			$regexp.="[ACGV]";
		}
		elsif ($_ eq "N"){
			$regexp.="[ACGUN]";
		}
		else{
			$regexp.=$_;
		}
	}
	return $regexp;
}


################################################## consolidation des fichiers CLADES et AFFILIATION
`cat $NGS_id_Results/$panam_output/PANAM_Affiliation_*.txt > $NGS_id_Results/$panam_output/PANAM_Affiliation.txt`;
`cat $NGS_id_Results/$panam_output/PANAM_Clades_NN_*.txt > $NGS_id_Results/$panam_output/PANAM_Clades_NN.txt`;
`cat $NGS_id_Results/$panam_output/PANAM_Clades_LCA_*.txt > $NGS_id_Results/$panam_output/PANAM_Clades_LCA.txt`;

#######################################################################################################


################################ reprise depuis 10-panam-clades.pl pour construction de %taxOTU
my %taxOTU;
my $c; my %hhh; my $bb;  my %bootstrap; my %sample; my %leveltax;

if (-e $NGS_id_Results."/".$panam_output."/PANAM_Affiliation.txt") { 
	open (AA, $NGS_id_Results."/".$panam_output."/PANAM_Affiliation.txt");

	
	while (<AA>) {	
		if ($_=~ /^>Read: (.*?)_/) {
			my $aa = $1;
			$sample{$aa}=1;
		}
		if ($_=~ /^>Read: (.*?)\t(.*?)\n/) {
			$c = $1; 
			$bootstrap{$c} = $2;
		}
		else {
			chomp($_);
			$hhh{$c}.= $_;
		}
	}
	close AA; 

	my %clade;my %nbr; my %tax; my %boot;

	foreach my $f (keys %hhh) { 
		if (($hhh{$f} !~ /Metazoa/) and ($hhh{$f} !~ /Other/)) {  ########### ajout du if le 9/10/2013
			$bb = $bootstrap{$f};
			my @tt = split (/Lowest node:/, $hhh{$f}) ;
			chomp ($tt[0]); chomp ($tt[1]);
		
		################"" %taxOTU pour user_file = pooled et %leveltax pour user_file = eachSample
	
			my @uu = split (/\t/, $tt[1]);
			my $node = $uu[0]; 

			$boot{'node'}{$node} = $bb;
			$clade{'node'}{$node}.=$f.", ";
			$nbr{'node'}{$node}++;
			$tax{'node'}{$node} =  $uu[1];
	
			if ((defined $uu[1])  and ($uu[1] ne "")) {
				$taxOTU{$f}{'LCA'} = $uu[1]; 

				$uu[1]=~s/\t//gi;
				$uu[1]=~s/^\s+//gi;
				$uu[1]=~s/\s+$//gi;

				my @aqs = split (/;/, $uu[1]);
				my $aqs = $#aqs; 

				#if ($aqs == 0) { 
				if ($aqs == 0) {}  ## modif 9/10/2013
					########" modif 23/5/2012
	# 						$taxOTU{$f}{'LCA'}=$uu[1]."Unclassified;";
					#$taxOTU{$f}{'LCA'}="Unclassified;";
					##########"fin modif
					#$leveltax{'LCA'}{'Unclassified'}++; 
				#}
				#elsif ($aqs == 1) {
				#	$taxOTU{$f}{'LCA'}=$uu[1]."unclassified ".$aqs[1].";";
				#	$leveltax{'LCA'}{$aqs[1]}{'unclassified '.$aqs[1]}++; 
				#}
				else {
				#	$taxOTU{$f}{'LCA'} = $uu[1];
					$leveltax{'LCA'}{$aqs[1]}{$aqs[2]}++;
				}
			}
	
			my @u = split (/\t/, $tt[0]);
			if ((defined $u[1])  and ($u[1] ne "")) {
				$taxOTU{$f}{'NN'} = $u[1];
				chomp($u[1]); 
			
				my @bqs = split (/;/, $u[1]); 
				my $bqs = $#bqs;

				#if ($bqs == 0) {
				if ($bqs == 0) {} # modif 9/10/2013
					###########"modif 23/5
	# 						$taxOTU{$f}{'NN'}=$u[1]."Unclassified;";
					#$taxOTU{$f}{'NN'}="Unclassified;" ;
					##############" fin modif
					#$leveltax{'NN'}{'Unclassified'}++; 
				#}		
				#elsif ($bqs == 1) {
				#	$taxOTU{$f}{'NN'}=$u[1]."unclassified ".$bqs[1].";";
				#	$leveltax{'NN'}{$bqs[1]}{'unclassified '.$bqs[1]}++;  
				#}
				else {
					$taxOTU{$f}{'NN'} = $u[1];
					$leveltax{'NN'}{$bqs[1]}{$bqs[2]}++
				}
			}

	###########################################"

			my @t = split (/Nearest neighbor:/, $tt[0]);
		
			my @proc = split (/\t/, $t[1]); 
			my $voisin = $proc[0]; 
	
			$boot{'voisin'}{$voisin} = $bb;
			$tax{'voisin'}{$voisin}=$proc[1]."\n\t\t";
			$clade{'voisin'}{$voisin}.=$f.", ";
			$nbr{'voisin'}{$voisin}++;
		} ### fin du if ajouté le 9/10/2013
	}

}
####################################### Fin reprise depuis 10-panam-clades.pl

my @dom_enlv = ("eukaryota", "bacteria" , "archaea" , "Metazoa", "Other");

foreach my $k (keys %taxOTU) { #print "k -- $k\n";
	foreach my $aff (keys %{$taxOTU{$k}}) { 
		foreach my $d (@dom_enlv) { 
			if (exists $dom{$d}) {}
			else {
				my @tab = split(/;/, $taxOTU{$k}{$aff}) ;
				#if ($taxOTU{$k}{$aff} =~ /$d/i) {  #print "$aff--- $k\t"; print "tax -- $taxOTU{$k}{$aff}\n";
				if ($tab[1] eq $d ) {
					delete ($taxOTU{$k}{$aff})
				}
			}
		}
	}
}
#################################
##### récupérer les OTUs avec seed, children et nombre ; selon le type d'analyse (Pooled ou eachSample)
##### récupérer $nbrSeq pour calcul des indices

my %nbrSeq;
my %seed; my %child;  
my %barIn;
my %NbrSeqClean;

if ($user_file eq "pooledSamples") { 
	open (FF, $NGS_id_Results."/".$preprocess_output."/pooled_sample/pooled_sample_OTU") || die "can not open file";
	
	#my %barIn;
	my %hps; 
	my $l;
	while (<FF>) {
		if ($_=~ /^OTU/) {}
		else {
			
			chomp ($_);
			my @tab = split (/\t/, $_);
			
			############################################ !!!!!!!!!!!
			if (exists $taxOTU{$tab[1]}{'LCA'}) {
				
				my $id = $tab[0]."\t".$tab[1];
				$hps{$id}{'id'} = $tab[1];
				my @tt = split (/, /, $tab[3]);			
				
				foreach my $e (@tt) {
				
					my $d;
					if ($e=~ /(\w.*?)_/) { 
						
						$d = $1; 
						########
					 	$NbrSeqClean{$d}++;
						#####
						$barIn{$d} = 1; 
						$hps{$id}{$d}++;
						$nbrSeq{$hps{$id}{'id'}}{$d}{'pooled'}++ ;
					}
				}	
		
				my $u = scalar(@tt);
				$u-= 1;
				my @t = @tt[0..$u];

				$child{$tab[1]} = \@t;
				foreach my $e (@t) {	
					if (defined $e) {
						$seed{$e} = $tab[1];
					}
				}
			} 
			else { delete $taxOTU{$tab[1]} }	### éliminer les seq qui n'ont pas été affiliées phylogénétiquement
		}
	}
	close FF;

	########################" OTU_distribution_taxo
	
	open (REPS, ">".$NGS_id_Results."/OTU_distribution_tax.txt");
 	print REPS "OTU_Id\tOTU_Seed\t";
	foreach my $pp ( (sort {$so{$a} <=> $so{$b}} keys %barIn)) { 									## VERSION 3.02 ##
		print REPS "$pp\t";
	}
	print REPS "NN taxonomy\tLCA taxonomy\n"; 
		
	foreach my $p (keys %hps) {
		print REPS "$p\t";
		foreach my $pp ( (sort {$so{$a} <=> $so{$b}} keys %barIn)) { 								## VERSION 3.02 ##

	 		if (defined $hps{$p}{$pp}) {
				print REPS"$hps{$p}{$pp}\t";
			}
			else {
				print REPS"0\t"
			}
		}
		print REPS "$taxOTU{$hps{$p}{'id'}}{'NN'}\t$taxOTU{$hps{$p}{'id'}}{'LCA'}\n" ;
	}
	close REPS;
}


if ($user_file eq "eachSample") {
	foreach my $s (keys %sample) { #print "s --- $s\n";
		open (FFF, $NGS_id_Results."/".$preprocess_output."/".$s."_OTU") || die "can not open file" ;
		while (<FFF>) {
			if ($_=~ /^OTU/) {}
			else {
				my @tab = split (/\t/, $_);
	
				############################################ !!!!!!!!!!!  /!\ NN ou LCA  dans le test d'existence ????
				if (exists $taxOTU{$tab[1]}{'LCA'}) {
					###### $nbrSeq
					$nbrSeq{$tab[1]}{$s}{'each'} = $tab[2]; #print "$s -- $tab[1] -- $tab[2]\n";
					#########			
				
					my @tt = split (/, /, $tab[3]);
					my $u = scalar(@tt);
					$u-= 1;
					my @t = @tt[0..$u];
					$child{$tab[1]} = \@t;
					foreach my $e (@t) {
						if (defined $e) {
							$seed{$e} = $tab[1]
						}
					}
				}else { delete $taxOTU{$tab[1]} }	### éliminer les seq qui n'ont pas été affiliées phylogénétiquement
			}
		}
		close FFF;
	}
}


##################################################################
	
my %seq_clean;					
foreach my $k (keys %taxOTU) { 
	foreach my $aff (keys %{$taxOTU{$k}}) { 	
		foreach my $e (@{$child{$k}}) {		
			if ((defined $e) and ($e ne "")) {
				my @oo = split (/_/,$e);			
				foreach my $s (%sample) { 										
					if ($oo[0] eq $s){ 										
						push (@{$seq_clean{$aff}{$s}}, $e) ;
					}
				}	
			}	
		}
	}
}

%child = (); 														## VERSION 3.01 ##

my %hash; my %min_norm;
my %min_norm_env;


$min_norm{'NN'} = $aff_norm;
$min_norm{'LCA'} = $aff_norm;

foreach my $aff (keys %seq_clean) {
	foreach my $ss(keys %{$seq_clean{$aff}}) {											# $ss = env 
		my $m = scalar (@{$seq_clean{$aff}{$ss}});
		if ($m >= $aff_norm) {
			$min_norm_env{$ss}=1;
		}
		foreach my $g (@{$seq_clean{$aff}{$ss}}) { 
			$hash{$ss}{$aff}{$g} = 1;
		}
	}
}

%seq_clean = ();													## VERSION 3.01 ##


##### sélectionner $min_norm séq dans chaque sample ####

my %subsample;
foreach my $k (keys %hash) { #print "k --- $k\n";										# $k = $env
	foreach my $aff (keys %{$hash{$k}}) { 
		my $sub=0;

		my $keys = keys(%{$hash{$k}{$aff}});		
		if ($keys >= $min_norm{$aff}){ #### éliminer les ech avec un nbr de seq < le seuil de normalisation
			foreach my $kk (keys %{$hash{$k}{$aff}}) { 
				$subsample{$k}{$aff}{$kk}=1;
				$sub++;
				if ($sub == $min_norm{$aff}) { last }
			}
		}
	}
}
#######

%hash = (); %min_norm = ();												## VERSION 3.01 ##		

my %d;
my %synt;
foreach my $k (keys %subsample) { 											# $k = $env
	foreach my $aff (keys %{$subsample{$k}}) {  
		foreach my $kk (keys %{$subsample{$k}{$aff}}) {
			my $aa = $seed{$kk};
			$synt{$k}{$aff}{$aa}{'reads'} .= $kk.",";
			$synt{$k}{$aff}{$aa}{'som'}++ ;
			$d{$aff}{$aa}= 1;
		}
	}
}	

%subsample = (); %seed = ();

my %distr; 
foreach my $v (keys %synt) {												# $v = $env
	foreach my $aff (keys %{$synt{$v}}) {
		foreach my $k (keys %{$d{$aff}} ) {
			if (defined $synt{$v}{$aff}{$k}{'som'}) {
				$distr{$k}{$aff}{$v} = $synt{$v}{$aff}{$k}{'som'}
			}
			else {
				$distr{$k}{$aff}{$v} = 0
			}
		}
	}
}

%d = (); 														## VERSION 3.01 ##		


if ($user_file eq "pooledSamples") { 
	open (ODNN, ">".$NGS_id_Results."/OTU_distribution_tax_normalized_NN.txt");
	open (ODNL, ">".$NGS_id_Results."/OTU_distribution_tax_normalized_LCA.txt");
	print ODNN "OTU_Seed\t"; print ODNL "OTU_Seed\t";
 
	foreach my $pp ( (sort {$so{$a} <=> $so{$b}} keys %synt)) { 								## VERSION 3.02 ##
		
		print ODNN "$pp\t" ;
		print ODNL "$pp\t" ;
	}
	print ODNN "NN taxonomy\n"; 
	print ODNL "LCA taxonomy\n"; 

	#### Print data
	foreach my $k (keys %distr) {
		print ODNL "$k\t"; print ODNN "$k\t";

		foreach my $p ( (sort {$so{$a} <=> $so{$b}} keys %synt)) { 								## VERSION 3.02 ##

			if (exists ($synt{$p})) {
			print ODNL "$distr{$k}{'LCA'}{$p}\t" ;
			print ODNN "$distr{$k}{'NN'}{$p}\t"
			}
			else {
			print ODNL "NA\t" ;
			print ODNN "NA\t"
			}
		}
		print ODNL "$taxOTU{$k}{'LCA'}\n";
		print ODNN "$taxOTU{$k}{'NN'}\n";
	}

	close ODNN;
	close ODNL;
}

%distr = ();  														## VERSION 3.01 ##		

	################################################ hashs pour calcul des indices et printing
	#########################################""

###################################################################################### BOUCLER sur les $env (%sample)
foreach $dbkey (keys %barIn) {
	print "$dbkey\n";
}
foreach my $env ( (sort {$so{$a} <=> $so{$b}} keys %barIn)) {
	
	chomp($env);															## VERSION 3.01 ##
	
	my %res;
	my %res_norm;															#if ($normaliz eq "yes") {
	#foreach my $env (keys %synt) {													# $v = $env
		foreach my $aff (%{$synt{$env}}) {
			foreach my $Id_OTU (%{$synt{$env}{$aff}}) {

				my @bqs = split (/;/, $taxOTU{$Id_OTU}{$aff});

				if ($#bqs == 0) {
					$res_norm{$aff}{$env}{$bqs[0]}{'nbrSEQ'}+= $synt{$env}{$aff}{$Id_OTU}{'som'};
					$res_norm{$aff}{$env}{$bqs[0]}{'OTU'}{$Id_OTU} = $synt{$env}{$aff}{$Id_OTU}{'som'}	;
				
					if ($synt{$env}{$aff}{$Id_OTU}{'som'} != 0) {
						$res_norm{$aff}{$env}{$bqs[0]}{'nbrOTU'}++;
					}
					if ($synt{$env}{$aff}{$Id_OTU}{'som'} == 1) {
						$res_norm{$aff}{$env}{$bqs[0]}{'1'}++;
					}
					if ($synt{$env}{$aff}{$Id_OTU}{'som'} == 2) {
						$res_norm{$aff}{$env}{$bqs[0]}{'2'}++;
					}
				}	

				else { 
					$res_norm{$aff}{$env}{$bqs[1]}{$bqs[2]}{'nbrSEQ'}+= $synt{$env}{$aff}{$Id_OTU}{'som'};
					$res_norm{$aff}{$env}{$bqs[1]}{$bqs[2]}{'OTU'}{$Id_OTU} =$synt{$env}{$aff}{$Id_OTU}{'som'};

					$res_norm{$aff}{$env}{$bqs[1]}{'nbrSEQ'}+= $synt{$env}{$aff}{$Id_OTU}{'som'};
					$res_norm{$aff}{$env}{$bqs[1]}{'OTU'}{$Id_OTU}= $synt{$env}{$aff}{$Id_OTU}{'som'};

					###
					if ($synt{$env}{$aff}{$Id_OTU}{'som'} != 0) {						
						$res_norm{$aff}{$env}{$bqs[1]}{$bqs[2]}{'nbrOTU'}++;
						$res_norm{$aff}{$env}{$bqs[1]}{'nbrOTU'}++;
					}

					if ($synt{$env}{$aff}{$Id_OTU}{'som'} == 1) {					
						$res_norm{$aff}{$env}{$bqs[1]}{$bqs[2]}{'1'}++;
						$res_norm{$aff}{$env}{$bqs[1]}{'1'}++				
					}
					if ($synt{$env}{$aff}{$Id_OTU}{'som'} == 2) {
						$res_norm{$aff}{$env}{$bqs[1]}{$bqs[2]}{'2'}++;
						$res_norm{$aff}{$env}{$bqs[1]}{'2'}++
					}
				}
			}
		}
	#}

########################### no normalizing

#if ($normaliz eq "no" ){
	# foreach my $env (keys %barIn) {																													# $env = $env
	foreach my $l (keys %taxOTU) { 
		my $nbrSeq; 
		if ($user_file eq "eachSample") {
			$nbrSeq = $nbrSeq{$l}{$env}{'each'};
		}
		if ($user_file eq "pooledSamples") { 
			$nbrSeq = $nbrSeq{$l}{$env}{'pooled'};
		}

		#print "$l --- nbrSeq --- $nbrSeq\n";
		################## parsing method ############## 

		foreach my $aff (keys %{$taxOTU{$l}} ) { 

			my @aqs = split (/;/, $taxOTU{$l}{$aff});
	
			if ($#aqs == 0 ) {
				$res{$aff}{$env}{$aqs[0]}{'nbrSEQ'}+= $nbrSeq;
				$res{$aff}{$env}{$aqs[0]}{'OTU'}{$l} = $nbrSeq	;
				if ($nbrSeq != 0) {
					$res{$aff}{$env}{$aqs[0]}{'nbrOTU'}++;
				}
				if ($nbrSeq == 1) {
					$res{$aff}{$env}{$aqs[0]}{'1'}++;
				}
				if ($nbrSeq == 2) {
					$res{$aff}{$env}{$aqs[0]}{'2'}++;
				}
			}	
			else { 
				$res{$aff}{$env}{$aqs[1]}{$aqs[2]}{'nbrSEQ'}+= $nbrSeq;
				$res{$aff}{$env}{$aqs[1]}{$aqs[2]}{'OTU'}{$l} = $nbrSeq;
				###
				$res{$aff}{$env}{$aqs[1]}{'nbrSEQ'}+= $nbrSeq;
				$res{$aff}{$env}{$aqs[1]}{'OTU'}{$l}= $nbrSeq;
				###
				if ($nbrSeq != 0) {
					$res{$aff}{$env}{$aqs[1]}{$aqs[2]}{'nbrOTU'}++;
					$res{$aff}{$env}{$aqs[1]}{'nbrOTU'}++;
				}
				if ($nbrSeq == 1) {
					$res{$aff}{$env}{$aqs[1]}{$aqs[2]}{'1'}++;
					$res{$aff}{$env}{$aqs[1]}{'1'}++;
				}
				if ($nbrSeq == 2) {
					$res{$aff}{$env}{$aqs[1]}{$aqs[2]}{'2'}++;
					$res{$aff}{$env}{$aqs[1]}{'2'}++;
				}
			}
		}
	}
	#}
#}

	######################################################### pas de normalisation ##################
	###################################### calcul 
	foreach my $parsing (keys %res) {
		foreach my $ss (keys %{$res{$parsing}}) { 
			foreach my $a (keys %{$res{$parsing}{$ss}}) { 
		
				my $nv2_n1; my $nv2_n2;

				################## calcul de schao
				
				$nv2_n1 = $res{$parsing}{$ss}{$a}{'1'};  $nv2_n2 = $res{$parsing}{$ss}{$a}{'2'};
			
				if ((($nv2_n1 >0) and ($nv2_n2 >= 0)) or (($nv2_n1==0 ) and ($nv2_n2 == 0))) {
					$res{$parsing}{$ss}{$a}{'chao'} =  $res{$parsing}{$ss}{$a}{'nbrOTU'} + (($nv2_n1*($nv2_n1 - 1))/(2*($nv2_n2+1)));
					$res{$parsing}{$ss}{$a}{'chao'} = sprintf("%.2f", $res{$parsing}{$ss}{$a}{'chao'});
				}
				else {
					$res{$parsing}{$ss}{$a}{'chao'} = $res{$parsing}{$ss}{$a}{'nbrOTU'} + (($nv2_n1*$nv2_n1)/(2*$nv2_n2));
					$res{$parsing}{$ss}{$a}{'chao'} = sprintf("%.2f", $res{$parsing}{$ss}{$a}{'chao'});
				}
		
				############### calcul de Coverage
				my $cov;
				if ($res{$parsing}{$ss}{$a}{'nbrSEQ'} != 0 ) {	
					$cov = 1- ($nv2_n1/$res{$parsing}{$ss}{$a}{'nbrSEQ'});
					$cov = $cov * 100;
					$res{$parsing}{$ss}{$a}{'cov'} =  sprintf("%.2f", $cov);
				}
				else {
					$res{$parsing}{$ss}{$a}{'cov'} =0
				}
		
			########################################"nv3

				foreach my $aa (keys %{$leveltax{$parsing}{$a}}) { 
					my $nv3_n1; my $nv3_n2;
		
				    ##################### calcul de schao
					$nv3_n1 = $res{$parsing}{$ss}{$a}{$aa}{'1'};  
					$nv3_n2 = $res{$parsing}{$ss}{$a}{$aa}{'2'};
									
					if ((($nv3_n1 >0) and ($nv3_n2 >= 0)) or (($nv3_n1==0 ) and ($nv3_n2 == 0))) {
						$res{$parsing}{$ss}{$a}{$aa}{'chao'} =  $res{$parsing}{$ss}{$a}{$aa}{'nbrOTU'} + (($nv3_n1*($nv3_n1 - 1))/(2*($nv3_n2+1)));
						$res{$parsing}{$ss}{$a}{$aa}{'chao'} = sprintf("%.2f", $res{$parsing}{$ss}{$a}{$aa}{'chao'});
					}
					else {
						$res{$parsing}{$ss}{$a}{$aa}{'chao'} = $res{$parsing}{$ss}{$a}{$aa}{'nbrOTU'} + (($nv3_n1*$nv3_n1)/(2*$nv3_n2));
						$res{$parsing}{$ss}{$a}{$aa}{'chao'} = sprintf("%.2f", $res{$parsing}{$ss}{$a}{$aa}{'chao'});
					}
			
				    ################## calcul de Coverage
					my $cov;
					if ($res{$parsing}{$ss}{$a}{$aa}{'nbrSEQ'} != 0 ) {	
						$cov = 1- ($nv3_n1/$res{$parsing}{$ss}{$a}{$aa}{'nbrSEQ'});
						$cov = $cov * 100;
						$res{$parsing}{$ss}{$a}{$aa}{'cov'} =  sprintf("%.2f", $cov);
					}
					else {
						$res{$parsing}{$ss}{$a}{$aa}{'cov'} =0;
					}
				}
			}
		}
	}

	###################################### printing
	foreach my $parsing (sort keys %leveltax) { 												# $parsing = LCA ou NN
		
		open (TAX, ">".$NGS_id_Results."/grp_tax".$parsing."_".$env."_tmp_all");
		print TAX "\n\t";
		
		foreach my $pp (sort keys %{$leveltax{$parsing}}) {
			if (($pp ne "") and ($pp ne ";")) {	
				print TAX "\n$pp\t";
				foreach my $ppp (sort keys %{$leveltax{$parsing}{$pp}}) {
					if (($ppp ne "") and ($ppp ne ";") and ($pp ne "Metazoa") and ($pp ne "Other")) {
						print TAX "\n\t$ppp";
					}
				}
			}
		}
		
		close TAX;

		my $o = $env;
				
		# foreach my $o (keys %{$res{$parsing}}) { 											## VERSION 3.01 ##
			#if ($o ne "") {
				open (OO, ">".$NGS_id_Results."/res_".$o."_".$parsing."_tmp_all");
				print OO "\t\t\t\t\t$o\n";
				print OO "\t#seq\t#OTUs\tSchao1\tShannon\tCoverage\n";
						
				foreach my $j (sort keys %{$leveltax{$parsing}}) { 
					if ((defined $j) and ($j ne "") and ($j ne ";")) {
						my $t2; my $sh=0; 
						if (exists $res{$parsing}{$o}{$j}{'nbrOTU'}) {
							foreach my $e (keys %{$res{$parsing}{$o}{$j}{'OTU'}}) {
								if ($res{$parsing}{$o}{$j}{'OTU'}{$e} != 0) {
									$t2=$res{$parsing}{$o}{$j}{'OTU'}{$e}/$res{$parsing}{$o}{$j}{'nbrSEQ'};
									$t2 = $t2 * log($t2);
									$sh += $t2;
								}
							}
				
							$sh = $sh*(-1);
							$res{$parsing}{$o}{$j}{'shannon'} = sprintf("%.2f", $sh);
							print OO "\t$res{$parsing}{$o}{$j}{'nbrSEQ'}\t$res{$parsing}{$o}{$j}{'nbrOTU'}\t$res{$parsing}{$o}{$j}{'chao'}\t$res{$parsing}{$o}{$j}{'shannon'}\t$res{$parsing}{$o}{$j}{'cov'}\n";
						}
							
						else {
							print OO "\t0\t0\t0\t0\t0\n";
						}
					}
					foreach my $jj (sort keys %{$leveltax{$parsing}{$j}}) { 
						my $sh=0; my $t2;
						if (exists $res{$parsing}{$o}{$j}{$jj}{'nbrOTU'}) {
							foreach my $e (keys %{$res{$parsing}{$o}{$j}{$jj}{'OTU'}}) {
								if ($res{$parsing}{$o}{$j}{$jj}{'OTU'}{$e} != 0) {
									$t2=$res{$parsing}{$o}{$j}{$jj}{'OTU'}{$e}/$res{$parsing}{$o}{$j}{$jj}{'nbrSEQ'};
									$t2 = $t2 * log($t2);
									$sh += $t2;
								}
							}
										
							$sh = $sh*(-1);
							$res{$parsing}{$o}{$j}{$jj}{'shannon'} = sprintf("%.2f", $sh);
		
							print OO "\t$res{$parsing}{$o}{$j}{$jj}{'nbrSEQ'}\t$res{$parsing}{$o}{$j}{$jj}{'nbrOTU'}\t$res{$parsing}{$o}{$j}{$jj}{'chao'}\t$res{$parsing}{$o}{$j}{$jj}{'shannon'}\t$res{$parsing}{$o}{$j}{$jj}{'cov'}\n";
						}
						else {
							print OO "\t0\t0\t0\t0\t0\n";
						}
					}	
				}
				close OO;
			#}
			`paste $NGS_id_Results"/grp_tax"$parsing"_"$o"_tmp_all" $NGS_id_Results"/res_"$o"_"$parsing"_tmp_all" > $NGS_id_Results/taxonomic_distribution_"$o"_"$parsing".txt`;

		#}																	## VERSION 3.01 ##

	## 		`paste  $NGS_id_Results"/"*"_tmp_all" > $NGS_id_Results/taxonomic_distribution_"$parsing".txt`;

		`rm $NGS_id_Results/*_tmp_*`;
	}
	
	%res = (); 																	## VERSION 3.01 ##
	 
	###############################################################


	######################################################### normalisation ##################
	###################################### calcul 
	foreach my $parsing (keys %res_norm) { 
		foreach my $ss (keys %{$res_norm{$parsing}}) { 
			foreach my $a (keys %{$res_norm{$parsing}{$ss}}) { 
		
				my $nv2_n1; my $nv2_n2;

				################## calcul de schao
				
				$nv2_n1 = $res_norm{$parsing}{$ss}{$a}{'1'};  $nv2_n2 = $res_norm{$parsing}{$ss}{$a}{'2'};
			
				if ((($nv2_n1 >0) and ($nv2_n2 >= 0)) or (($nv2_n1==0 ) and ($nv2_n2 == 0))) {
					$res_norm{$parsing}{$ss}{$a}{'chao'} =  $res_norm{$parsing}{$ss}{$a}{'nbrOTU'} + (($nv2_n1*($nv2_n1 - 1))/(2*($nv2_n2+1)));
					$res_norm{$parsing}{$ss}{$a}{'chao'} = sprintf("%.2f", $res_norm{$parsing}{$ss}{$a}{'chao'});
				}
				else {
					$res_norm{$parsing}{$ss}{$a}{'chao'} = $res_norm{$parsing}{$ss}{$a}{'nbrOTU'} + (($nv2_n1*$nv2_n1)/(2*$nv2_n2));
					$res_norm{$parsing}{$ss}{$a}{'chao'} = sprintf("%.2f", $res_norm{$parsing}{$ss}{$a}{'chao'});
				}
		
				############### calcul de Coverage
				my $cov;
				if ($res_norm{$parsing}{$ss}{$a}{'nbrSEQ'} != 0 ) {	
					$cov = 1- ($nv2_n1/$res_norm{$parsing}{$ss}{$a}{'nbrSEQ'});
					$cov = $cov * 100;
					$res_norm{$parsing}{$ss}{$a}{'cov'} =  sprintf("%.2f", $cov);
				}
				else {
					$res_norm{$parsing}{$ss}{$a}{'cov'} =0
				}
		
			########################################"nv3

				foreach my $aa (keys %{$leveltax{$parsing}{$a}}) { 
					my $nv3_n1; my $nv3_n2;
		
				    ##################### calcul de schao
					$nv3_n1 = $res_norm{$parsing}{$ss}{$a}{$aa}{'1'};  
					$nv3_n2 = $res_norm{$parsing}{$ss}{$a}{$aa}{'2'};
									
					if ((($nv3_n1 >0) and ($nv3_n2 >= 0)) or (($nv3_n1==0 ) and ($nv3_n2 == 0))) {
						$res_norm{$parsing}{$ss}{$a}{$aa}{'chao'} =  $res_norm{$parsing}{$ss}{$a}{$aa}{'nbrOTU'} + (($nv3_n1*($nv3_n1 - 1))/(2*($nv3_n2+1)));
						$res_norm{$parsing}{$ss}{$a}{$aa}{'chao'} = sprintf("%.2f", $res_norm{$parsing}{$ss}{$a}{$aa}{'chao'});
					}
					else {
						$res_norm{$parsing}{$ss}{$a}{$aa}{'chao'} = $res_norm{$parsing}{$ss}{$a}{$aa}{'nbrOTU'} + (($nv3_n1*$nv3_n1)/(2*$nv3_n2));
						$res_norm{$parsing}{$ss}{$a}{$aa}{'chao'} = sprintf("%.2f", $res_norm{$parsing}{$ss}{$a}{$aa}{'chao'});
					}
			
				    ################## calcul de Coverage
					my $cov;
					if ($res_norm{$parsing}{$ss}{$a}{$aa}{'nbrSEQ'} != 0 ) {	
						$cov = 1- ($nv3_n1/$res_norm{$parsing}{$ss}{$a}{$aa}{'nbrSEQ'});
						$cov = $cov * 100;
						$res_norm{$parsing}{$ss}{$a}{$aa}{'cov'} =  sprintf("%.2f", $cov);
					}
					else {
						$res_norm{$parsing}{$ss}{$a}{$aa}{'cov'} =0
					}
				}
			}
		}
	}

	###################################### printing
	foreach my $parsing (sort keys %leveltax) {
		open (TAX, ">".$NGS_id_Results."/grp_tax".$parsing."_".$env."_norm_tmp_all");
		print TAX"\n\t";
		
		foreach my $pp (sort keys %{$leveltax{$parsing}}) {
			if (($pp ne "") and ($pp ne ";")) {	
				print TAX "\n$pp\t";
				foreach my $ppp (sort keys %{$leveltax{$parsing}{$pp}}) {
					if (($ppp ne "") and ($ppp ne ";") and ($pp ne "Metazoa") and ($pp ne "Other")) {
						print TAX "\n\t$ppp";
					}
				}
			}
		}
		
		close TAX;
		
		#foreach my $o (keys %{$res_norm{$parsing}}) {
		#foreach my $o (keys %barIn) { 
		my $o = $env;
			#if ($o ne "") {	
				open (OO, ">".$NGS_id_Results."/res_".$env."_".$parsing."_norm_tmp_all");
				print OO "\t\t\t\t\t$o\n";
				print OO "\t#seq\t#OTUs\tSchao1\tShannon\tCoverage\n";
						
				foreach my $j (sort keys %{$leveltax{$parsing}}) { 
					if ((defined $j) and ($j ne "") and ($j ne ";")) {
						my $t2; my $sh=0; 
						######
						if (exists $res_norm{$parsing}{$o}) {
						#####
							if (exists $res_norm{$parsing}{$o}{$j}{'nbrOTU'}) {
								foreach my $e (keys %{$res_norm{$parsing}{$o}{$j}{'OTU'}}) {
									if ($res_norm{$parsing}{$o}{$j}{'OTU'}{$e} != 0) {
										$t2=$res_norm{$parsing}{$o}{$j}{'OTU'}{$e}/$res_norm{$parsing}{$o}{$j}{'nbrSEQ'};
										$t2 = $t2 * log($t2);
										$sh += $t2;
									}
								}
				
								$sh = $sh*(-1);
								$res_norm{$parsing}{$o}{$j}{'shannon'} = sprintf("%.2f", $sh);
								print OO "\t$res_norm{$parsing}{$o}{$j}{'nbrSEQ'}\t$res_norm{$parsing}{$o}{$j}{'nbrOTU'}\t$res_norm{$parsing}{$o}{$j}{'chao'}\t$res_norm{$parsing}{$o}{$j}{'shannon'}\t$res_norm{$parsing}{$o}{$j}{'cov'}\n";
							}
							
							else {
								print OO "\t0\t0\t0\t0\t0\n";
							}
						}
						else {
							print OO "\tNA\tNA\tNA\tNA\tNA\n";
						}


					######
					}
					foreach my $jj (sort keys %{$leveltax{$parsing}{$j}}) { 
						my $sh=0; my $t2;
						######
						if (exists $res_norm{$parsing}{$o}) {
						#####
							if (exists $res_norm{$parsing}{$o}{$j}{$jj}{'nbrOTU'}) {
								foreach my $e (keys %{$res_norm{$parsing}{$o}{$j}{$jj}{'OTU'}}) {
									if ($res_norm{$parsing}{$o}{$j}{$jj}{'OTU'}{$e} != 0) {
										$t2=$res_norm{$parsing}{$o}{$j}{$jj}{'OTU'}{$e}/$res_norm{$parsing}{$o}{$j}{$jj}{'nbrSEQ'};
										$t2 = $t2 * log($t2);
										$sh += $t2;
									}
								}
										
								$sh = $sh*(-1);
								$res_norm{$parsing}{$o}{$j}{$jj}{'shannon'} = sprintf("%.2f", $sh);
		
								print OO "\t$res_norm{$parsing}{$o}{$j}{$jj}{'nbrSEQ'}\t$res_norm{$parsing}{$o}{$j}{$jj}{'nbrOTU'}\t$res_norm{$parsing}{$o}{$j}{$jj}{'chao'}\t$res_norm{$parsing}{$o}{$j}{$jj}{'shannon'}\t$res_norm{$parsing}{$o}{$j}{$jj}{'cov'}\n";
							}
							else {
								print OO "\t0\t0\t0\t0\t0\n";
							}
						}
						else {
							print OO "\tNA\tNA\tNA\tNA\tNA\n";
						}
					}	
				}
			#}
			close OO;
			`paste $NGS_id_Results"/grp_tax"$parsing"_"$o"_norm_tmp_all" $NGS_id_Results"/res_"$o"_"$parsing"_norm_tmp_all" > $NGS_id_Results/taxonomic_distribution_"$o"_"$parsing"_norm.txt`;

		#}

	## 		`paste  $NGS_id_Results"/"*"_tmp_all" > $NGS_id_Results/taxonomic_distribution_"$parsing".txt`;

		`rm $NGS_id_Results/*_tmp_*`;
		
	}
	 
	###############################################################

	foreach my $k (keys %NbrSeqClean) {
		print "$k --- $NbrSeqClean{$k}\n"
	}

}

my @normtypes = ("", "_norm");
my @parsingmeth = ("LCA", "NN");

my $head = `head -n 1 $NGS_id_Results/OTU_distribution_tax.txt`;
my @samples = split("\t", $head);
shift(@samples); shift(@samples);
pop(@samples); pop(@samples); 
chomp(@samples);

foreach my $parsing (@parsingmeth) {
	foreach my $ntype (@normtypes) {

		my $cpt = 1;
		
		foreach my $env (@samples) {	
			my $coden = "";
			if ($cpt<10) {
				$coden = "0";
			} 
			#print $env."\n";
			if (!( -e "tax$parsing$ntype"."_tmp0" ) )  {
				`awk -F\"\t\" '{print \$1"\t"\$2;}' $NGS_id_Results/taxonomic_distribution_$env"_"$parsing$ntype.txt > $NGS_id_Results/tax$parsing$ntype"_"tmp0`;
			}
			`awk -F\"\t\" '{print \$3"\t"\$4"\t"\$5"\t"\$6"\t"\$7"\t"\$8;}' $NGS_id_Results/taxonomic_distribution_$env"_"$parsing$ntype.txt > $NGS_id_Results/tax$parsing$ntype"_tmp"$coden$cpt`;
			`rm $NGS_id_Results/taxonomic_distribution_$env"_"$parsing$ntype.txt`; 
			$cpt++;
		}
		`paste $NGS_id_Results/tax$parsing$ntype"_"tmp* > $NGS_id_Results/taxonomic_distribution_$parsing$ntype.txt`;
		`rm $NGS_id_Results/tax$parsing$ntype"_"tmp*`;
	}
}


